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^ ; Abstract 
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■ We have simulated the classical Heisenberg antiferromagnet on a triangular 

os : 

■ lattice using a local Monte Carlo algorithm. The behavior of the correlation 

ctf ; 

^ ■ length £, the susceptibility at the ordering wavevector x(Q)> an d the spin 

-A ■ 

stiffness p clearly reflects the existence of two temperature regimes - a high 

o 

O ■ temperature regime T > T^, in which the disordering effect of vortices is 

> : 

dominant, and a low temperature regime T < T^, where correlations are con- 
trolled by small amplitude spin fluctuations. As has previously been shown, 
in the last regime, the behavior of the above quantities agrees well with the 
predictions of a renormalization group treatment of the appropriate nonlinear 
sigma model. For T > T t h, a satisfactory fit of the data is achieved, if the tem- 
perature dependence of £ and x(Q) i s assumed to be of the form predicted 
by the Kosterlitz-Thouless theory. Surprisingly, the crossover between the 
two regimes appears to happen in a very narrow temperature interval around 
T th ~ 0.28. 
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I. INTRODUCTION 



Magnetic ordering phenomena of both classical and quantum antiferromagnets on non- 
bipartite lattices are a fascinating subject. The simplest and most frequently studied model 
of this type is the Heisenberg antiferromagnet on a triangular lattice (HAFT). While the 
question of whether or not the typical features of this model have been observed in exper- 
iments is still a controversial issue M, the theoretical understanding of these features has 
advanced rapidly during recent years In contrast to antiferromagnets on bipartite 

lattices, the HAFT exhibits non-collinear magnetic order in its classical and most likely 
also in its quantum ground state. As a consequence, the order parameter of the HAFT is 
represented locally by a set of three mutually orthogonal unit vectors or, alternatively, by a 
rotation matrix which defines the local orientation of this set relative to some fixed frame of 
reference. Renormalization group (RG) studies of appropriate nonlinear sigma (NLer) mod- 
els H|J have revealed a number of interesting properties of the HAFT. The symmetry of the 
model was found to be dynamically enhanced from 0(3) (g> 0(2) to 0(4), and in a two loop 
RG calculation for the classical HAFT ||, the temperature dependence of the correlation 
length £ was obtained as 



£ = A Ci G yT/B e B ' T , (1) 

where A is the lattice constant, B = \/3ir(y + |) = 6.994. The prefactor C^ RG is left 
undetermined by the RG calculation. 

It follows from topological considerations that the order-parameter field of the HAFT 
allows for excitations of the form of Z 2 vortices |7]f|. A numerical study of the classical 
HAFT has revealed that these vortices become abundant above a threshold temperature 
T t h ~ 0.3 and that they unbind for T > T t h, similarly as the Z vortices in the planar 
XY model above the Kosterlitz-Thouless (KT) transition temperature T^t 0- Further 
evidence for this similarity beween the dissociation mechanism of the Z 2 vortices and that 
of the Z vortices has been provided by Kawamura and Kikuchi |T(J. In recent work ||TT 



we 



have studied the influence of the vortices on the partition function of the classical HAFT 
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on the basis of the NLcr model. While a true KT-type phase-transition can be ruled out 
for the HAFT, our results suggest that for T > T t h the vortices will affect the properties 
of the HAFT rather drastically. In particular, the disorder induced by the unbinding of 
the vortices can be expected to lead to a crossover from the T-dependence Eq. ([]]) of the 
correlation length in the low temperature regime T < Tth to a KT-type behavior 

£ = AC* KT eMb/(T-T th ) l 2) (2) 

in the high temperature regime. 

It is the aim of the present paper to supplement our recent analytical study ||11|| , which 



was based on a continuum description of the HAFT by a Monte Carlo (MC) simulation of 
the original lattice model. A similar study has recently been published by Southern and 



Young |12|| . While we closely follow the method of these authors, our conclusions will be 
quite different from theirs. 

In the next Section, we first give a brief account of the technique we used. Subsequently 
we present the results for the correlation length £ and the antiferromagnetic susceptibility 
x(Q). While these quantities are directly accessible to simulations in the high temperature 
regime, where the disordering effect of the vortices limits the range of the correlations, the 



key quantity to be computed in the low temperature regime is the spin stiffness ||13| , |T2 
Our numerical results for this last quantity will be presented and discussed in Section III. 
Finally, we summarize the evidence for a vortex induced crossover transition in Section IV. 



II. CORRELATION LENGTH AND ANTIFERROMAGNETIC SUSCEPTIBILITY 

The classsical Hamiltonian of the triangular Heisenberg antiferromagnet can be defined 

as 

ff=ESrS 3 . (3) 

<ij> 

Here, the Sj are three dimensional unit vectors and the sum extends over all distinct pairs of 
nearest neighbor sites of a triangular lattice of L 2 sites. The exchange constant has been set 
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to unity. The classical ground state of the Hamiltonian Eq. is a coplanar arrangement 
in which the spins on the three sublattices are oriented at 120° relative to each other, 

Si = eicos(QRj) + e 2 sm(QRj) . (4) 

Here, e 1; e 2 are a pair of mutually orthogonal unit vectors and Q can be any one of the 
six vectors pointing towards the corners of the hexagonal Brillouin zone of the triangular 
lattice, e.g. Q = ^(|,0). The correlation length £ can be obtained assuming a Ornstein 
Zernicke form for the structure factor 

5(q) = -^E eiq(R ^ Rj) < S ^ S , > ( 5 ) 

S{Q) 

i+e(q-Q) 2 ' 

x(Q) = S(Q) /T is then the susceptibility of the system at the ordering wavevector. The spin 
correlations < Sj ■ Sj > in Eq. (pp can be determined by MC techniques. We used the local 
algorithm described by Kawamura and Miyashita J7|. The lattice is divided into independent 
sublattices. Then, the spins of each of these sublattices are updated sequentially. For a given 
spin, a new direction is chosen at random and the standard Metropolis rule is used to decide 
whether the new direction is to be accepted. If it is discarded, a precessional motion through 
a randomly chosen angle about the direction of the local mean field is performed. We apply 
this method to systems of linear sizes L — 12 • 2 n ,n = 0, 1, ..5. For the smaller systems, 
n < 3, we discard the first 2 • 10 4 sweeps for equilibration and average over the next 2 • 10 5 
sweeps. For n = 4, we average over 4 • 10 5 sweeps after discarding the initial 10 5 sweeps, 
and for n = 5, the average is over 1.8 • 10 6 , and 2 • 10 5 sweeps are discarded. 

A selection of results for the correlation length £ and for the antiferromagnetic suscep- 
tibility x(Q) which have been obtained by averaging over 3-5 independent runs of these 
lengths is tabulated in Table |. As will become apparent shortly, the data shown in this 
Table are crucial in checks of theoretical predictions for the temperature dependence of £ 
and x(Q). To exhibit possible finite size effects, the Table contains two pairs of data for 
each temperature which correspond to two different system sizes L, 2L. In general, £ and 
x(Q) decrease with T, but increase with the system size L. If, for a given temperature 
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T and a given system size L , the data for £ and x(Q) exhibit no size dependence upon 
doubling the system size, then one can conclude that the system size L Q suffices to obtain 
size independent data for all T > Tq. Data which are size independent by this criterion are 
marked by an asterisk in Table Q. Obviously, we cannot exclude that the data for the lowest 
temperature T = 0.3 obtained for the L = 384 system are still size dependent. Certainly, 
however, our data for T = 0.3 are lower bounds to the thermodynamic limits of £ and x(Q) 
at this temperature. To facilitate the comparison of our data with the RG predictions we 
include in Table | the values for £ and x(Q) which result from fits of the expressions Eq. 
(H|) and Eq. (§) to these data @. 

In Figs. H] and |2|, we show our complete sets of results for the correlation length and 
for the antiferromagnetic susceptibility as functions of the temperature. Obviously, for any 
given system size L, there is an inflection point in the sequences of data for £ and x(Q). 
This point defines a temperature Tl below which both £ and x(Q) begin to exhibit finite 
size effects. In fact, as can be seen in Fig. [I], the correlation length increases linearly with 
the system size for sufficiently low temperatures T <C T L . Figs. [I] and |2] also contain fits 
of different theoretical predictions to the numerical data. The dashed lines represent fits of 
the RG result, Eq. ([]]), to our data for £ and of the form 

X (Q)=C* G (T/B) 4 eM2B/T), (6) 



proposed by Southern and Young |L2| on the basis of RG calculations, to our data for 
x(Q). In these RG predictions, the constants C RG and C\ G are the only undetermined 
parameters. In our fits, we neglect the data points for temperatures T < Tl which contain 
finite size effects. In agreement with Southern and Young [12 we find C RG ~ 3 • 10 -8 and 



C\ G ~ 6 • 10 12 . The solid lines represent fits of the KT forms Eq. (0) and 

x(Q) = C x KT exp( 7 -.b/(T-T th ) l i) (7) 

to the data. With the KT form for £, Eq. (^), the last expression follows from the general 
relation S(Q) ~ £ 2_?? for the structure factor at the ordering wave vector, if r\ is assumed to 
take the value rj = 1/4 as for a proper KT transition. For the threshold temperature, we use 



the value T t h = 0.28, which we can be inferred from the temperature dependence of the spin 
stiffness as will be discussed in the next section. This leaves the constant b which is common 
to the expressions Eq. (fj) and Eq. ([F]) and the constant factors C^ KT and C\ T of Eq. (0) 
and Eq. (|7|) as fit parameters. As before, we ignore data points for T < Tl in our fits. The 
results are b = 0.77, C% T = 0.47 and C\ T = 0.40. Obviously, for temperatures T > 0.3, the 
exp(6/(T — T t / l )5)-temperature dependence of the KT forms fits the data better than the 
exp(_B/T)-temperature dependence predicted by the RG analysis. In Figs. § and [| we plot 
£ and x(Q) logarithmically against (T — T^) -1 / 2 so that the KT forms Eq. (fj) and Eq. (^) 
appear as straight lines. These lines are seen to fit the data quite well in the temperature 
interval 0.30 < T < 0.34, whereas the agreement between the curves representing the RG 
forms is restricted to a narrow interval around T « 0.31. In particular, we emphasize that 
for T = 0.3, the RG predictions are incompatible with the data points for £ and x(Q) which 
are listed in Table | with their respective errors. In this context, we recall that if our data 
for T = 0.3 do not represent the thermodynamic limits of £ and x(Q); they are certainly 
lower bounds to these limits. Hence, the discrepancy between the RG predictions and the 
true values of £ and x(Q) may even be larger than has been inferred here. We note that 
the fit of the KT form, Eq. (J?]), to the data for the susceptibility x(Q) is better than that 
of the RG form, Eq. @, to £. This may be attributable to the lower quality of the data for 
£ which are obtained indirectly from the Ornstein-Zernicke expression Eq. (^) in the limit 

eiQ-q| «i- 

In summary, we observe that our results combined with the earlier findings of Kawamura 
and Miyashita support the view that in the temperature range T > T t h the spin corre- 
lations of the HAFT are decisively influenced by unbound vortices, so that a perturbative 
treatment of the model is inadequate in this temperature regime. It should be obvious, how- 
ever, that in the above analysis of our numerical results we have been guided by our previous 



prediction |]TTJ] that, in the case of the HAFT, the vortex unbinding mechnism leads to a 
KT-type temperature dependence of the correlation length above a crossover temperature 
T t h- While we do not claim to have found compelling evidence for this prediction, we regard 
our numerical results as strong support for it. 
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III. SPIN STIFFNESS 



To further corroborate the above view and in order to get insight into the low temperature 
regime, where the correlation length exceeds the accessible system sizes, we also determined 
the spin stiffness p in our simulations. 

The diagonal components p ai a = 1,2,3, of the spin stiffness tensor are the second 
derivatives of the free energy density f(9 a ) with respect to the twist angles 8 a of the spins 
around three mutually orthogonal axes e a |i3|Jl2]| , 



2 

<i,j> 



< 



^ ] Sj x Sj ■ e Q (u ■ e, 

<i,j> 



2 



> • (8) 



VSTL 2 

Here, u is the lattice direction along which the twist is applied and is the direction of 
the bond between nearest neighbor lattice sites i and j. The prefactor in Eq. (H) has been 
chosen such that Eq. (Jj|) is the stiffness per unit area. 

In the simulation, the thermal averages on the right hand side of Eq. (K) are replaced 
by averages over configurations which are generated by a large number of successive MC 
sweeps. For a finite system, the ordered spin structure will change its orientation in space 
in the course of the simulation. In a sufficiently large number of sweeps, one will therefore 
measure the average spin stiffness 

P =\iip a - (9) 

6 a=l 

Since there is no long range order in the HAFT for any finite temperature, p must vanish 
for the HAFT in the thermodynamic limit for any finite temperature. However, for finite 
system sizes L, p will be finite for sufficiently low temperatures such that £ > L. According 
to Eq. (|I|), this condition should be satisfied for system sizes L < 10 3 up to temperatures 
of the order of unity, unless the constant C^ RG is exceedingly small as has been suggested by 
Southern and Young W%. As we have argued above, the results for £ in the high temperature 



regime, T > T t h, should not be interpreted as evidence for such a small value of C^ RG . In 
fact, a naive integration of the 2-loop renormalization group equations which starts with 



the microscopic parameters of the HAFT as initial values yields the result Eq. ([!]) with 

In our simulations, we determine the three p a separately in each sweep and thus obtain 
the averages p for each sweep. In Fig. [5], we show the average spin stiffness as a function of 
the temperature for system sizes L = 12, 24, . . . 384. The steep drop in p which occurs as T 
increases beyond T t h — 0.28 is consistent with the rapid decrease of £ in the same temperature 
regime in which the vortices become unbound [0. In contrast to the spin stiffness of the 
planar XY model, p does not saturate in the low temperature regime T < T t h with increasing 
system size L but decreases with increasing L. This behavior is to be expected, since in 
contrast with the correlation length of the XY model, the correlation length of the HAFT 
remains finite for low temperatures, where the vortices are bound in pairs. As T — > 0, our 
data for p approach the correct limiting value 1 / v^3- 

In Fig. H we display p for various temperatures and system sizes. The data are averages 
over 3-5 independent runs of lengths comparable to those which we have described above, 
error bars indicate the standard deviation. 

In the low temperature regime, the thermodynamics of the classical HAFT should be 
captured by the appropriate NLa model OHffl]. On the basis of a renormalization group 



treatment of this model, Azaria et al. |1| have made detailed predictions for the dependence 
of the spin stiffness on the linear system size L and the correlation length £. In their MC 



study of the classical HAFT, Southern and Young |L2| found excellent agreement with the 
predicted L dependence of the spin stiffness tensor at the temperature T = 0.2. 

In order to be able to compare our numerical results with the predictions of the RG 



analysis of the NLcr model, we integrated the two loop RG equations of Azaria et al. [{15 
starting from the initial conditions p(L = A) = I/a/3 and ps(L = A)/px(L = A) = 2. 
Here, p\ and p3 are the two main components of the spin stiffness tensor with respect to 



the reference frame of the local order parameter [0. By the above initial conditions we 
identify p\ and p3 with their microscopic values on the scale of the lattice constant A. We 
find that in the temperature regime under consideration, T < 0.3, the average stiffness p 
varies linearly with In L to a very good approximation on the scale 12 < L < 384 covered by 
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our simulations, see Fig. [7| below. From the RG equations, one can also infer that the slope 
p'{L) = —dp{L)/d\aL decreases from p' = T/(3tt) to p' = T/(4ir), when L increases from a 
value of the order of the lattice constant, L ~ A, to a value of the order of the correlation 
length, L ~ £. In Fig. ^ the straight lines are least squares fits to the data. The slopes of 
these lines, normalized to the maximal theoretical slope p' = T / (37r), are tabulated in Table 
0. For T < 0.25, the slopes are seen to be rather close to their maximal value which obtains, 
when L is of the order of the lattice constant. This is not unexpected since according to 
the RG calculations, the correlation length is many orders of magnitude larger than our 
maximal system size L = 384 for these temperatures. The larger values of p' which we find 
for T > 0.28 are incompatible with the RG theory. Hence we conclude that for T > 0.28 
unbound vortices, not being taken into account by the RG analysis, begin to limit the range 
of the spin correlations in the HAFT. 

If one defines the correlation length £ through the matching condition p(L = £) = 0, 
then the result of the integration of the two loop RG equations can be cast into the following 
form 

p = f{T,L) ln(f/L) . (10) 

This is shown for two different temperatures in the two graphs in Fig. |7]. Obviously, the 
function f(T, L) depends weakly on L. The relation Eq. (p~0| ) makes it possible to obtain 
the correlation length from the Monte Carlo simulation, even in the low temperature regime 
where £ is much larger than the system size L. Inserting our MC data for p into Eq. ( |T0"D 
and solving for £, we obtain the data points shown in Fig. |8] for T < 0.28. This Figure also 
includes the data for T > 0.28 which have already been shown in Fig. |l|. The dashed and 
solid lines in Fig. | represent fits of the RG and KT forms, Eqs. (1) and (2), to the MC 
data for T < 0.28 and T > 0.28, respectively. 

IV. SUMMARY 

We have shown that the MC simulation of the HAFT yields compelling evidence for 
the influence of vortices on the spin correlations and hence on the thermodynamics of this 
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model. In agreement with earlier findings by Kawamura and Miyashita |7j we find that the 
disordering effect of the vortices sets in rather abruptly at a temperature T t ^ ~ 0.28. Up to 
this temperature, our data for the spin stiffness p and the ensuing temperature dependence 
of the correlation length £ are consistent with the predictions of the RG analysis of the 
HAFT which ignores the existence of topological defects such as vortices. For T > T t h, 
however, the simulation yields temperature dependences of the correlation length £ and 
the antiferromagnetic susceptibility x(Q) which are incompatible with the RG predictions. 
Instead, in this temperature regime, the temperature dependences of £ and x(Q) which 
follow from the KT picture of unbinding vortex pairs provide satisfactory fits of the data. 
A rapid increase of the density of unbound vortices for T > 0.3 had already been found 
by Kawamura and Miyashita in their simulation of the HAFT 0. It had not been clear, 
however, whether this phenomenon would lead to the same temperature dependences of the 
correlations of the HAFT as had been predicted for the planar XY model by Kosterlitz and 
Thouless 0. In our recent analytical study [PH| , we were led to the conclusion that this 
should indeed be the case. The present numerical study fully supports this conclusion. As 



we have discussed in Ref. [IT|, the crossover transition from the RG type behavior to the 
KT type behavior of the correlations of the HAFT cannot imply a phase transition in the 
proper sense, because the correlation length is finite both below and above the transition 
temperature T^. The results shown in Fig. |8] indicate, however, that the transition happens 
in a very narrow interval around T^. Therefore we consider it possible that the derivative 
of the correlation length with respect to the temperature develops a discontinuity at T t ^ in 
the thermodynamic limit. 
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TABLES 

TABLE I. Correlation length £ and antiferromagnetic susceptibility x(Q) for system sizes 
L = 96, 192, and 384 for various temperatures. The statistical error represents the standard 
deviation over 3-5 independent runs. The two last columns contain the RG predictions for £ and 



x(Q). Asterisks indicate agreement of the results for different system sizes, see main text. 



T L 


£ X 


£,RG XRG 


0.300 384 
0.300 192 


105.0 ± 11.9 5014 ± 543 
83.5 ± 3.5 2800 ± 135 


oZ.o oDUy 


0.305 384 
0.305 192 


60.7 ± 6.7 2184 ± 273 

* 

63.4 ± 3.8 1899 ± 55 


O/.U 1/3D 


0.310 192 
0.310 96 


41.3 ± 1.8 1040 ± 39 
38.2 ± 0.3 748 ± 20 


on v 01 A 


0.315 192 
0.315 96 


31.7 ± 1.1 594 ± 25 
29.3 ± 0.7 494 ± 11 


28.0 476 


0.320 192 
0.320 96 


22.7 ± 2.0 349 ± 27 

* 

22.6 ± 1.1 350 ± 13 


19.9 253 


0.330 192 
0.330 96 


14.3 ± 3.4 169 ± 17 

* 

14.5 ± 1.5 169 ± 9 


10.4 76 


0.340 192 
0.340 96 


11.8 ± 3.3 103 ± 7 

* 

9.8 ± 1.1 96 ± 5 


5.7 25 


TABLE II. Slopes of the fits in Fig. |?j normalized to the maximum slope T/(3ir). 


T 

P 1 Pmax 


0.10 0.20 0.25 0.28 0.29 0.30 
0.962 0.937 0.939 1.045 1.191 2.098 
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FIGURES 

FIG. 1. Correlation length £ as a function of the temperature T. Solid line: the KT form 
Eq. d) with C% T = 0.47, b = 0.77, and T th = 0.28. Dashed line: the RG behavior Eq. © with 
4g = 3-10- 8 - 

FIG. 2. Antiferromagnetic susceptibility x(Q) as a function of the temperature T. Solid line: 
the KT form Eq. (0) with C£ T = 0.40, 6 = 0.77, and = 0.28. Dashed line: the RG behavior 
Eq. © with C* G = 6-10- 12 . 

FIG. 3. Correlation length £ as a function of (T — T^)"2 . The lines are defined as in Fig. [|. 

FIG. 4. Antiferromagnetic susceptibility x(Q) as a function of (T — T t ^) - 2. The lines are 
defined as in Fig. |2[ 

FIG. 5. Spin stiffness p as a function of the temperature T. 

FIG. 6. Spin stiffness p as a function of the system size L. The straight lines are mean squares 
fits, their slopes are tabulated in Tab. [n|. 

FIG. 7. Spin stiffness p as a function of ln(£/L) as obtained from the two loop RG equations. 
The curves correspond to T = 0.1 and 0.25, respectively. 

FIG. 8. Correlation length as a function of the temperature. Dashed lines: fit of the RG result 
to the low temperature data. Solid line: fit of the KT form to the high temperature data. The 
inset shows the data on a larger temperature scale. 
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